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This article deals with a viscoplastic material model of overstress type. The model 
is based on a multiplicative decomposition of the deformation gradient into elastic 
and inelastic part. An additional multiplicative decomposition of inelastic part is 
used to describe a nonlinear kinematic hardening of Armstrong-Frederick type. 

Two implicit time-stepping methods are adopted for numerical integration of 
evolution equations, such that the plastic incompressibility constraint is exactly 
satisfied. The first method is based on the tensor exponential. The second method 
is a modified Euler-Backward method. Special numerical tests show that both ap- 
proaches yield similar results even for finite inelastic increments. 

The basic features of the material response, predicted by the material model, are 
illustrated with a series of numerical simulations. 
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Notation 

F deformation gradient 

Fi inelastic part of the deformation gradient 

Fe elastic part of the deformation gradient 
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Fii dissipative part of Fi 

Fie energy storage part of Fi 

/C current configuration 

K, reference configuration 

/C stress-free intermediate configuration 

IC intermediate configuration of microstructure 

C right Caucliy-Green tensor (see (j6])) 

Ci inelastic right Cauchy-Green tensor (see ([7]) J 

Cii inelastic right Cauchy-Green tensor of microstructure (see 

Ce elastic right Cauchy-Green tensor (see 

Cie elastic right Cauchy-Green tensor of microstructure (see ([8]) 2) 

E Green strain tensor (see ([9])^) 

F Almansi strain tensor (see 

T Cauchy stress tensor 

S weighted Cauchy tensor (Kirchhoff tensor) (see (fT9l) ) 

S, T 2nd Piola-Kirchhoff tensors operating on /C, /C, respectively (see ( |20|) 1 

X, X, X backstress tensors operating on /C, /C and respectively (see (12^ ) 

H the driving force for inelastic flow (see fl37|) ^) 

S the driving force for inelastic flow of microstructure (see fl37j) n) 

1 second-rank identity tensor 

M* , M* covariant pull-back and push- forward (see (jl]) ) 
(M~^)*, (M^^)* contravariant pull-back and push-forward (see ^) 

A ■ B = AB product (composition) of two second-rank tensors 

A : B scalar product of two second-rank tensors (see (!22|) ) 

II A II I2 norm of a second-rank tensor (Frobenius norm) (see (HOl) j^) 

II A II* induced norm of a second-rank tensor (spectral norm) (see (1701) ) 

(■)° deviatoric part of a tensor (see (HOl) o) 
transposition of a tensor 
inverse of transposed 

tr(-) trace of a second-order tensor 



A 
(■) 



covariant Oldroyd rate with respect to /C (see (IT^ J 



covariant Oldroyd rate with respect to K, (see (fT5|) n) 

(■) unimodular part of a tensor (see 

sym(-) symmetric part of a tensor (see ([H 

skew(-) skew-symmetric part of a tensor (see fl67|l o) 

{x) MacCauley bracket (see (HT]) n) 

ip specific free energy 

5i specific internal dissipation (see fl3Tl) ) 

K initial yield stress 

R isotropic hardening 

*i? trial isotropic hardening (see 

s inelastic arc length 

Sd dissipative part of s 
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Se 


energy storage part of s (see (|T^) 


Ai 


proportionality factor (inelastic multiplier) (see (j^Tl)]^) 


/ 


overstress (see fj^T|)n) 


Sym 


space of symmetric second-rank tensors 




norm of the driving force (see ([S3D) 




incremental inelastic parameter (see 0781)) 




mass density in the reference configuration 


k 


bulk modulus (see (126|)) 




shear modulus (see fl26l)) 


c 


bulk modulus of microstructure (see (123) i) 
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hardening modulus (see (|27l)^) 




technical strain (see (1931)^^) 




shear strain (see flMl)]^) 




axial and shear stresses, respectively 



1 Introduction 



New materials, such as ultrafine-grained- aluminium (see the papers [15], [2H]), are of 
special interest for many practical applications. To promote the innovation of the new 
materials, the robust numerical simulation of the material response is required. It is de- 
sirable to have a phenomenological description of the material which on the one hand 
takes important phenomena into account, and on the other hand enables stable numerical 
computations. 

In this paper we investigate the simulation of rate-dependent material behavior with 
equilibrium hysteresis effect (for the general introduction to the theory of viscoplasticity 
see, for example, [33], [23], [T2]). 

The Bauschinger effect is observed in most metals under non-monotonic loading. The 
most popular approach to describe the Bauschinger effect was proposed by Armstrong 
and Frederick pj in 1966. Application of the Armstrong-Frederick hardening concept 
within the framework of Perzyna type viscoplasticity (see [33]) yields the classical material 
model of overstress type (see [3], [23]). This model has the advantage that it admits 
simple rheological interpretation (see fig. [T]a). Such phenomena as creep, relaxation and 
nonlinear kinematic hardening are taken into account by the model. Simple modification 
of this model is possible to include isotropic hardening as well HI . 

Several strategies can be adopted for the generalization of this model to finite strains (see. 



^ The diagram in fig. [TJa provides insight into the rheological modeling of kinematic hardening. 
To the best of our knowledge, there is no simple rheological diagram of viscoplastic material 
with isotropic hardening. 
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for example, [6], [38], [26], [37], |2l], [13], [31]). Some of the generalizations were analyzed 
numerically in [5]. Following the elegant approach of Lion [2l], we use the rheological 
interpretation (fig. [TJa) of the classical model to construct its finite-strain counterpart. 

The specific assumptions of the material modeling used in this paper are as follows: 

• Multiplicative decomposition of the deformation gradient into elastic and inelastic part: 
F = F,Fi ([21], [22]). 

• Multiplicative decomposition of the inelastic part into energy storage part and dissipa- 
tive part: F; = FieFji ([23). 

• Free energy is a sum of appropriate isotropic strain energy functions ([24]). 

The resulting material model takes both kinematic and isotropic hardening into account. 
The thermodynamic consistency is proved. 

The purpose of the present paper is threefold. First, we formulate the material model 
under consideration. In particular, we transform the constitutive equations to the reference 
configuration in order to simplify the numerical treatment. Next, two implicit schemes 
for the numerical integration of evolution equations are developed. Finally, we analyse 
numerically the basic properties of the material response, predicted by the model. 

A global implicit time stepping procedure in the context of displacement based FEM re- 
quires a proper stress algorithm (local integration algorithm) [40]. Such algorithm provides 
the stresses and the consistent tangent operator as a function of the strain history locally 
at each integration point. A set of internal variables is used in this paper to describe the 
history dependence, and the stress algorithm includes implicit integration of a system of 
differential (evolution) and algebraic equations. 

Two most popular implicit schemes for integration of inelastic strains in the context of 
viscoplasticity/plasticity are: 

• Backward- Euler scheme, also referred as implicit Euler scheme (see, for example, |10j . 
|35J, [36J, il3j, ^). 

• Exponential scheme, also referred as Euler scheme with exponential map (see, for ex- 
ample, [3S], [21], [30], 0). 

The exponential scheme is advantageous since it retains the inelastic incompressibility 
even for finite time steps. Thus, an important geometric property of the solution is au- 
tomatically preserved. Moreover, the numerical error of Euler-Backward method, related 
to the violation of incompressibility, tends to accumulate over time (see, for example, [5], 
[T4]). Therefore, even for small time steps, the numerical solution deviates from the exact 
solution after some period of time. 

Helm [13] modified the classical Euler-Backward scheme, using a projection on the group 
of unimodular tensors, to enforce the incompressibility of inelastic flow. 
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In this work we implement in a uniform manner both modified Euler-Backward method 
(MEBM) and the exponential method (EM). Both methods result in a nonlinear system 
of equations with respect to strain-like internal variables Ci = Fi'"Fi, Cii = FjfFii and 
^ = Ai At 0. This nonlinear system is split into two subproblems: 
First subproblem: Finding Ci, Cii with a given ^. 

Second subproblem: Finding ^, such that an incremental consistency condition is satisfied. 
This adapted strategy is more robust than the straightforward application of a nonlinear 
solver to the original system of equations. At the same time, this approach is not limited by 
the special form of the free energy, and finite elastic strains are likewise allowed. Moreover, 
the stress algorithms are applicable in the limiting case of rate-independent plasticity (as 
viscosity tends to zero). 

Although the material response is anisotropic, it is shown that MEBM as well as EM 
exactly preserve the symmetry of Ci and Cii. Furthermore, the accuracy of both integration 
algorithms is verified with the help of special numerical tests. Both methods provide 
similar results with almost the same integration error. A common feature of MEBM and 
EM is that the numerical error is not accumulated over time. 

The phenomenological description of each specific material can be schematically subdi- 
vided into three steps: 

• Material testing, such that the important phenomena make themselves evident. 

• Choosing an appropriate phenomenological model, that reproduces qualitatively the 
experimental data. 

• Parameter identification, using the experimental data. 

To illustrate the basic characteristics of the material model we simulate a series of material 
testing experiments. These experiments are uniaxial tension and torsion under monotonic 
and cyclic loading. In particular, we conclude that the material model can be used (after 
a proper parameter identification) to describe the mechanical response of an aluminium 
alloy processed by ECA-pressing [15], [28] . 

Throughout this article, bold-faced symbols denote first- and second-rank tensors in M^. 
Expression a := b means a is defined to be another name for b. 



2 Material model of finite viscoplasticity 

The material model is motivated by the rheological diagram in fig. [Ha. This diagram 
takes the kinematic hardening of Armstrong-Frederick type into account (for the sake of 
simplicity the isotropic hardening is omitted in the diagram). The total inelastic strains 
and the inelastic strains of microstructure are used as internal variables. The evolution of 
these quantities is closely related to the energy dissipation during the inelastic processes. 

^ ^ > is an incremental inelastic parameter, defined by (|78p 
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Besides, additional real-valued strain-like internal variables are introduced in order to 
describe a nonlinear isotropic hardening. 



2.1 Kinematics 



For a fixed time instant t > let /C C be a current configuration occupied by the solid. 
Suppose /C C is the reference configuration, which uniquely designates the material 
points. Let us consider the motion law in the form : /C ^ /C. For every point P G /C 

we define the deformation gradient tensor F := ' . The deformation gradient F 

transforms a material line element ciX on the reference configuration /C into a current 
material line element (ix 

dx = F dX. (1) 

Let us consider the classical multiplicative decomposition of the deformation gradient F 
into elastic part Fg and inelastic part F; ( |21j 



F = FeFi. (2) 

The mechanical justification uses the idea of the local (within a neighborhood of the mate- 
rial point) elastic unloading. The transformation rule ([T]) is represented as a combination 
of two linear operators 

c/x= Fe(Fi rfX). 

Therefore, we can interpret Fj ciX as a fictitious material line element on some intermedi- 
ate configuration fC (see fig.lUb). We will call this configuration the stress-free intermediate 
configuration. 

A second multiplicative decomposition is introduced in order to simulate a nonlinear 
kinematic hardening of Armstrong- Frederick type. Following Lion [21], we decompose the 
inelastic part Fj into energy storage part Fie and dissipative part Fa 

Fi = FieFii. (3) 

The energy storage part Fje describes the heterogeneity of elastic strains associated with 
the energy storage on the microscale. The dissipative part Fji can be attributed to slip 
processes on the microscale (see [23], [13] for details). Decomposition ([3]) implements the 
intermediate configuration of micro structur^^ K, (see fig. [Hb). The commutative diagram 
in fig. [Hb summarizes both multiplicative decompositions. 

In this paper we deal with second-order tensors, which operate on configurations /C, /C, /C, 
/C. Pull-back (push-forward) operations describe the transformation of tensor fields during 
the change of configurations. Let M G |f, F;, Fii, Fe, Fjej be a linear transformation of 



In [24j the similar configuration is called the intermediate configuration of kinematic hardening. 
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energy 
' dissipation ' 



energy 
storage 



microstructure 



energy 
dissipation 



energy 
storage 




Fig. 1. Rheological model (a) and a commutative diagram (b) showing corresponding configu- 
rations with transformations of material line elements. 

material line elements on two different configurations. We define corresponding pull-back 
and push-forward of covariant tensor field by 



M*(-) := M^(-)M, M,(-) := M-^(-)M-\ 
Pull-hack and push-forward of contravariant tensor field are given by 
(M-^)*(-) = M-i(-)M-^, (M-T),(-) = M(-)M^. 
Thus, the right Cauchy-Green tensor is a covariant pull-back of 1 

C := F*l = F^F. 



(4) 



(5) 



(6) 



In the same manner, we define the inelastic right Cauchy-Green tensor Ci and inelastic 
right Cauchy-Green tensor Cii of microstructure: 



Q := F*l = F^Fi, Cii := F^l = FjFii. 



(7) 



Further, the elastic right Cauchy-Green tensor Ce and the elastic right Cauchy-Green 
tensor Cie of microstructure are defined by 



Ce := F:1 



C = F* 1 = F F 

^le • ■'-10-'- -■- le-"- le- 



The tensors 

E:=^(C-1), r := F,E = F-^EF-^ (9) 

are termed the Green strain tensor and the Almansi strain tensor, respectively. Basing on 
E, we define a corresponding strain tensor F, which operates on fC 

f := (Fi),E = FrTEFri. 

Multiplicative decomposition ([2]) implements the additive decomposition of F: 

F = FeFi f = fi + fe, (10) 



7 



where f i is a purely inelastic Almansi tensor 

1 



and Fe is the elastic Green tensor 

fe:=^(FjFe-l) = ^(Ce-l). (H) 

Analogously, multiplicative decomposition ([3]) implements the additive decomposition of 
the pull-back of the inelastic Almansi tensor Fi to JC: 

f ■ •= F* f ■ = F^f -F- 

■•- 1 • ■•- le-*- 1 ■•- le-*- 1-"- 165 
Fi = FieFii fi = fii + fie, (12) 

where 

f , := ^(1 - FrrjFr,'), t-,, := ^{fIF-,, - 1) = i(Qe - 1). (13) 
Finally, we define the inelastic pull-back of F; to /C 

r^ := F*f i. 

In this paper, the evolution of isotropic hardening is taken into account, similar to the 
Armstrong-Frederick rule. To this end, we introduce two real-valued internal variables 
of strain type: s and s^- The first variable is the classical inelastic arc length, and Sd is 
interpreted as a dissipative part of s, such that 

Se := s - Sd (14) 

controls the energy stored due to the isotropic hardening (see section 2.3). 



2.2 The concept of dual variables 



The formalism of dual variables developed by Haupt and Tsakmakis [TT] specifies the 
choice of stress and strain variables as well as their time derivatives. According to this 

A 

concept, we introduce the covariant Oldroyd rates (■), (■) with respect to the stress-free 
configuration /C and the microstructural configuration JC, respectively, 

A:=Fi.((F*(-))'), (0:=Fii.((Fr,(-))'), (15) 
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where ()' stands for material time derivative. The alternative representation of covariant 
Oldroyd derivatives is as follows 

0)=(-)- + Lf(-) + (-)Li, (0=(-)- + L^(-) + (-)Lii, (16) 
Li := F,Fr\ L, := F,,Fr:\ 

Equation (fT5|) yields the inelastic deformation rates Fi and Ta as symmetric parts of Lj 
and Lii, respectively: 

A 

f i= sym(Li), f ii= sym(Lii), (17) 
where the symmetric part of a tensor is given by 

sym(-) := ^ 



Let T be the Cauchy stress tensor. The weighted Cauchy tensor (or Kirchhoff stress 
tensor) is defined by 

S := (detF)T. (19) 

Now, we define the 2nd Piola-Kirchhoff tensors operating on K. and /C using a contravariant 
pull-back of S 

S:=(F;T)*S, t:=(F-T)*S=(Fr^)*S. (20) 

The introduced stress and strain tensors form the following conjugate pairs: (S, F), (S, F), 
(T, E), such that the work and the stress power are invariant under the change of config- 
uration: 

A 

S : F = S : f = t : E, S :f = t : E. (21) 
Here " : " denotes the scalar product of two second-rank tensors 

A: B :=tr(A-B^). (22) 



Next, we denote by X the backstress tensor, which operates on the intermediate configura- 
tion )C of microstructure. This tensor can be interpreted as a generalized force, associated 

<> 

with strain measure Fi and strain rate Fi. According to the concept of dual variables, 
we define transformation rules for the backstress tensor, such that quantities X : Fj and 

X : Fi remain invariant under the change of configuration: 

X := (Fi:T)*X = F,1xF,t, X := (Fi-^)^X = F,,±FI (23) 

A 

X : f i = X : f i = X : fi, X : f i= X :f i= X : f;. (24) 
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2.3 Free energy 



Suppose that the free energy is given as a sum of isotropic function£3 (cf. [21], [13]) 

= V^(f e, fie, Se) = ^e\{T,) + V^kin(f ic) + V'iso(Se), (25) 

where the tensors fe and fie are defined by (fTTl) . (fT3l) o. Here, '?/'ei(f e) corresponds to the 
energy, stored due to macroscopic elastic deformations. The "inelastic" part ^/'kin(fie) + 
'?/'iso(se) represents the energy, stored in the microstructure during the viscoplastic flow 
due to the heterogeneity of dislocations. The following special form of the free energy can 
be used 

PnM^.) = ^(ln\/detCe)' + ^(trS^- 3), (26) 

PR^kin(rie) = ^(trCic - 3j, PRViso(Se) = " {Sef ■ (27) 

Here, fc>0,/i>0,c>0,7>0 are material constants. The overline (■) denotes the 
unimodular part of a tensor 

A := (det A)-^/3A . (28) 

Denote by ^'^^^ the derivative of real-valued function a with respect to tensor-valued 
argument A such that 

Sa=^4^: SA. (29) 

Using this notation, we introduce formally the following potential relations for stresses S, 
X and for isotropic hardening R 

^ dlpelite) ^ (9V'kin(fie) p 5^iso(Se) ,„nA 

S = p.^j— , X = p,,^p^, i? = PK.^^. (30) 



2.4 Clausius-Duhem inequality 



The Clausius-Duhem inequality imposes an additional constraint on the material response, 
which states that the internal dissipation is always nonnegative. For isothermal processes 
the specific internal dissipation 5i takes the form (see |T2] ) 

(5i := — t : E-^/- > 0. (31) 

Now, let us rewrite this expression, using relations of previous subsections. First, we note 
that S and X are isotropic functions of Fg and fie, respectively. In particular, since S and 

^ the first two terms of this additive split can be motivated by the rheological model (fig. [TJa). 



10 



Fe commute, we get 

S : (Li^fe + feLi) # 2(f eS) : 2(f ,S) :f ^ . (32) 

Further, note that 

S-t^S :f e +S :f S :f e +S : {Ijt, + t^i) + S 1;^ S :f e +(CeS) 1; . (33) 
In the same way, since X and Fie commute, we obtain from (|T2l) . (fT6i) n. (fT7l) n 

X:fi=X:fie+(CieX) :fii. (34) 
Thus, we get for the stress power 



S :F ^ S :F -X :Fi +X :Fi 



A . A 



S :Fe +(CeS - Xj :Fi +X :Fie +(CieXj :Fii . (35) 
Hence, the internal dissipation takes the form 

5, = : E - ^ Is :f # (U - ^) :f . +(^X 

pR Pr Vr 5Fe^ 

A 

:f 

Pr ' ' pR 

We abbreviate 



+ -(CeS - X) :f i +-(QeX) :f , ie. (36) 



C^V'kin 


■) :f ie 


:fii- 









S := CeS - X, H := CieX. (37) 

Finally, taking into account potential relations fl30l) and definition f[T^ . we simplify fl36l) 
to obtain the Clausius-Duhem inequality in the form 

pJi = S :Fi -i? s + H :Fii + i? > 0. (38) 



2.5 Evolution equations 



Following the standard procedure, we formulate the evolution equations for internal vari- 
ables so that inequality fl38l) holds for arbitrary mechanical loadings (cf. [23], [T3]). 

Ti = Ai--^, f ii:= Ai X H°, (39) 
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-Ai, Sd := -sR, 
3 7 

1 



||A|| := VA : A, A° := A - -tr(A)l, (40) 

3 

where the inelastic multipher Ai is determined according to the Perzyna rule [32] 

1 / 1 \ ™ - [2 r 1 

Ai:=-(;^/) , /:=||S°||-^-[ir + i?], (x):=max(x,0). (41) 

Here, x>0,r]>0,m>l, P>0, K>0 are material parameters, fco > is used to get 
a dimensionless term in the bracket. 

Let us show that inequality (1381) is fulfilled. For instance, we prove that the parenthetical 
term in (138|) is nonnegative. Indeed, 



sli -R ^)=Ai(||S°l|-w|i?) 



JO if / < (Ai = 0) 
Ai(/ + ./fir) if / > (Ai > 0) 



Therefore, the material model, defined in sections 2.3 and 2.5, is thermodynamically con- 
sistent. 

According to the evolution equations fl39l) . the tensors and S° are termed the driving 
force for inelastic flow and the driving force for inelastic flow of micro structure, respec- 
tively. Note that both fiows are incompressible. In fact, 

(detFi)' = (detFi) tr(^ f i ^ = 0, (detFii)' = (detFii) tr(^ Fii ^ = 0. (42) 

Under appropriate initial conditions, it follows from fH2]) that 

detFi = detFii = detCi = detCii = 1. (43) 

The inelastic flow takes place if the overstress / is positive. A case of rate-independent 
plasticity is covered by these evolution equations as viscosity 77 tends to zero. 



2.6 Transformation to the reference configuration 



A direct numerical treatment of Oldroyd derivatives in fl39|) . formulated with respect 
to flctitious configurations /C and is complicated. In this subsection we rewrite the 
material model in terms of strain-like internal variables Ci, Cii, s, s^ such that the rate 
of Ci, Cii will be given by the material time derivatives Ci, Cii. The transformation of 
the model includes: 

• representation of the free energy %Ij through C, Ci, Cii, s, s^. 
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transformation of the potential relations for stresses, 
representation of through C, Cj, Cii. 

transformation of the evolution equations. 



2.6.1 Representation of the free energy 

Let (Ji, J2, J3) be a full system of invariants of a second-rank tensor, defined by 
Ji(A) := tr A, MA) := ^tr A^, J,(A) := hi Al 

Then, using multiplicative decompositions ([2]), ([3]) and the property 
tr(AB) = tr(BA), it is easily proved that 

</fe(Ce) = Jfc(CCi "^), Jfc(Cie) = Jfc(CiCii '^), k = 1,2,3. 

Since V'ei and ipkin are isotropic functions, it follows from (l25l) . (144|) that 

^jj = V^(C, Ci, Cii, S, Sd) = TZ-ellCCi^^) + V^km(CiCii-l) + ?/^i,o(s - Sd). 

S.^.S Transformation of the potential relations for stresses 
Recall that (cf. equation (9.60) in [T2]) 



T^a(ABA^ _ 9a(ABAT) 



A=const 



a(ABAT) dB 
On the other hand. 

But, 

Substituting F;"'^, C, Ce and F^J , Ci, Cje for A, B, ABA"^ in (gS]), we obtain 



~_ gv^ei(cCi-^) ~_ dijUCiC^r'] 

J- — -^Pr „ ^, ^ — ^Pr- 



dC 



Ci=const (^Ci 



Cii=const 



Now let us show that CT and CiX are isotropic functions of CCi ^ and CjCji ^, 
tively. Indeed, note that 

daiAB-') _^«(AB_^3-i^ a4(A) ^ fc = l,2,3. 



dA 



B=const 



diAB-') ' dA 
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Combining (H6|) and (H7|) . we get 



^1 dJk{CC, 



T:TT(CiCii 



If the special form (1261) . (127|) for the free energy is used, then, taking into account 
pressibihty relations (H3|) . we obtain 



T = A; InJdet(C) C"^ + C-^(CCi 



Representation of WE^W 
First, let us note that 

trS = tr(CeS - X) = tr(Cf - QX). 
Next, we compute the inelastic contravariant and covariant pull-back of S° 



= Cr^CT - X - trS Cr^ ^ Cri(ct - QX 



~ X D 



(Fi)*S° = CtCi - CiXCi - trS Q ^ (ct - CiX)°Ci. 

Furthermore, since S° G Sym, we get 

llSDf = tr(S°S^)=tr|[(FrT)*sD] [(Fi)*S°]}. 
Substituting fl50l) and fl5T!) in fl52!) we obtain the norm of the driving force 



^ = ^tr[(ct- CiX 



2.6.4 Transformation of the evolution equations 
Note that 

tr H = tr (CicX) = tr(CiX) 
Next, we compute the covariant pull-back of : 



Covariant Pull-back of ( l39ll yields 



Q = 2(fO* m 2^(fO C. = 2(F,)* t-^ 2A.x(f,)*H- 



1(39) 1 
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Table 1 

Summary of the material model 



Ci = 2|(CT-CiX)°Ci, 



Cii = 2Aix(CiX)^Cii, 



Ci|i=o = CP, detCO = l, 
Cnk=o = COi, detCOi = l, 



s : 



1 - Ipn — ^ — 



s\t=Q = S°, Sd|i=0 = Sd, 



I Ci=const ' 
R = JSe, Se = S - Sd, 



I Cii=const ' 



tr 



(CT - CiX) 



Combining this with 0511) and (1551) . we obtain 



Q = 2|(ct - QX)°Ci, Qi = 2Aix(QX)°Cii. 



(56) 



The material model is summarized in table [H 



3 Integration algorithms 

The exact solution of fl56|) has under proper initial conditions the following geometric 
property: Ci, Cii he on the manifold M, defined by 

M := {b G 5?/m : detB = l}. (57) 

Hence, system fl56p is a system of differential equations on the manifold (cf. the paper 
[H]). In this section we analyse two numerical schemes, such that the numerical solution 
lies exactly on M. 

3.1 Modified Euler- Backward and exponential scheme 

Consider the Cauchy problem for a system of nonlinear ordinary differential equations 

A{t) = f(A(t),t)A(t), A(0) = A°, det(A°) = 1. 
Suppose that the tensor-valued function f is sufficiently smooth, tr(f(B,t)) = and 

(f(B,t))^B G Sym yBeSym, k = 1,2,3,.. . . (58) 
Under such conditions the exact solution lies on M. 
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Remark: condition (l58l) is nontrivial, since f(B,t) is, in general, an anisotropic function 
of B. 

By "A,"'^^A denote numerical solutions respectively at t„ and tn+i. 

At := — t„. Suppose that "A G M is given. The classical Euler-Backward method 

(EBM) uses the equation with respect to the unknown """"^A : 

A =[l-At f ("+^A, tn+i)] "A. (59) 

Recall that for small B 

[1-B]~^ = 1 + B + B2 + B3 + ... . (60) 
The exponential method (EM) is based on the equation 

"+iA = exp (At f ("+^A, t„+i)) "A, (61) 
where the tensor exponential is given by 

exp (b) := 1 + B + ^B^ + Ib^ + ... . (62) 

Let us show that both methods yield a symmetric solution. The idea of the proof is as fol- 



lows. Substituting dHUD for [l-At f ("+iA, t„+i)J in (^,and(^ for exp i^At f C'+^A, t„+i 
in (l6T]) . we get for both methods 

"+^A = K("+^A), (63) 



K(B):= l + Atf(B,Wi) + ^Cfc(f(B,t„+i)) "A, (64) 

^ k=2 ^ 

with some coefficients c^. Next, let us consider an auxiliary problem 

Aaux = sym(K(Aaux))- (65) 

Here, the symmetrization operator sym(-) = \{[-) + (O""^^ is used. Suppose Aaux is a 
solution of (|65|) . According to properties (|58|) . since Aaux is symmetric, we obtain 

Aaux "A-^ Aaux € Sym, K(Aaux) Aaux e Sym. (66) 

Subtracting (l66ll ^ from ( lUBll o and taking (1651) into account, we get 

skew(K(Aaux)) "A-i Aaux G Sym, skew(-) := (■) - sym(-). (67) 
Here skew(-) stands for the skew-symmetric part of a tensor. Thus, (1671) ^ yields 

skewfskew(K(Aaux)) "A-^ Aaux) = 0. (68) 



The Neumann series (j60p converges if ||B||* < 1. 
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Since skew skew(-) = skew(-), from (168|) follows 



skew(K(Aaux)) = skew^skew(K(Aauxy 
Denote by || ■ ||* an induced norm of a tensor 



1 - " A-i A 



|B| 



max ||Bx||2, \\'x\\2 -J xl + + x^. 

||xj|2 = l 



Then, 

||AB||* < ||A||*||B||* 
Note also that for small At 

111 -"A-i A 



||skew(A)||* < ||A| 



(69) 

(70) 

(71) 
(72) 



Taking the norm of both sides of (l69l) and using flTTl) , fl72|) , we get 

l|69] l. l(7T|l o 



||skew(K(Aa 



< 



skew(K(Aa 



1 - "A-i A 



< ||skew(K(Aaux))r ||l-"A-i Aauxir < -||skew(K(Aa 



This implies that ||skew(^K(Aaux)) ||* = 0. Therefore, K(Aaux) ^ Sym and Aaux is a 
solution of (!63|) . In other words, equations (!63l) and (l65l) are equivalent. ■ 

This means that no modifications of (!59|) anc? (!6T!) are necessary to ensure the symmetry 
of the solution "+^A. The reader will have no difficulty in showing that the problem of 
symmetry does not occur also for a system of equations of type (IMj) . (1641) . 



The following modifications of equations (1591) and (I6T1) leave the corresponding original 
solutions unchanged: 



n+l 



A = sym{ [l - At f ("+i A, t„+i)] ' "A}, 
"+iA = sym{ exp (At f ("+1A, t„+i)) "A}. 



(73) 
(74) 



The advantage of the exponential method based on (1611) or (1741) is that the constraint 
det("'"''^A) = 1 is exactly satisfied. In this paper we modify the right-hand side of fl59l) 
and (I73|) . using the projection (■) = (det(-))~^/^(-) on the group of unimodular tensors 
(cf. [14]): 



1-At f("+iA,t„+i) 



^A, 



(75) 



A = sym{ [l - At f ("+i A, t„+i)] ^ "A}. (76) 

Let us remark that both (1751) and (!76|) yield the same solution ""*"^A G M. Thus, the 
modified Euler-Backward (MEBM) is formulated by (1751) or (1761) . Further, we notice that 



17 



the exponential method (EM) (l6Ti) is equivalent to 



n+1 



A = sym{ exp (At f ("+iA, t„+i)) "A}. 



(77) 



Remark. We have a freedom in choosing between fl75|) and fl76|) for MEBM. Similarly, 
the EM can be based either on fl6T|) or fl77j) . In this paper we use symmetrized equations 
dZSD, frrr|) . The reason is that these two equations can be formulated with respect to six 
real unknowns. At the same time equations (1751) . fl^T]) are formulated with respect to nine 
independent real unknowns. 

3.2 Adaptation of integration methods to the evolution equations 

Suppose that the deformation gradient "+^F at the time t^+i = tn + At is known. Further, 
assume that the internal variables Ci, Cii, s, Sd at the time t„ are given by "Ci, ""Cji, "s, "sa, 
respectively. In this subsection we formulate a system of equations for finding the internal 
variables at the time tn+i- 

First, we adopt the modified Euler-Backward scheme (1761) and the exponential scheme (1771) 
to the numerical integration of evolution equations (1561) . We stress that the right-hand 
sides in ( 15B1) satisfy requirements (155]) . For instance, let us analyse the evolution equation 
for Ci. Note that, since CT and CjX are isotropic functions of CCi^^ and CiCii"^, 



The evolution equations for s, Sd are discretized by implicit Euler scheme. Further, con- 
sider an incremental inelastic parameter 



2^(ct - QX) = rfil + d^CCr^ + c/3(CCr')' + d^C.Cr,^ + 4(CiC 
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-1\2 



with some suitable dn G M. It remains to check that 



(CCr^)'=Ci G Sym, (Q 



Cr^)''Ci e Sym, V C;, Qi G Sym, k = 1,2,3,... . 



^ := At "+^Ai. 



(78) 



Finally, we get the following system of equations. 



n+1 



Ci - sym(Ki("+iC,"+iCi,"+iCii,0) = 0, 



(79) 




(80) 




(81) 
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n+1 n 



/5 2, 



= y tr 

where the operators K^, k G {i, ii} are defined by 

^1 



~ \ D' 



1-B, 



exp 



"Cfc if MEBM is employed 
*Cfc if EM is employed 



~ X D 



~ \ D 



(82) 
(83) 
(84) 



^5) 
^6) 



Bii("+iCi,"+^Cii,0 := 2 e x("+^Ci "+'X) 

Here "+it, "^^X are functions of "+^0, "+^Ci, "+^0;;, given by (06]) (or by (gH]) if the 
special form of ipei, ipkin is used). 



3.3 Solution strategy 



First, we exclude "+is,"+^Sd from ([82]), (I83])i (cf. [Hj) to get 



^+^R = R{0:= ^^y^^ ^i^-TC^s-'^Sd). (87) 



Next, substituting ([HID and ([HTj) in ([HSDg, we represent ^^+7 as a function of ^^+^0;, "+^Cii, ^. 
Thus, the problem is reduced to system (1791) . (180|) . (18T!) with respect to """"^Ci, "+^Cii, ^. 

In this paper we decompose problem (!79|) . (IHOj) . (IHTl) as follows. The variables ""'"^Ci, ""'"^Cii 
are uniquely determined by system (!79|) . (IHO!) with a given ^. Let us denote the correspond- 
ing solution by (qC'+^C, 0, Cii("+^C, O) • Substituting this solution in ([H]), we obtain 
a function ^{''+^0,^). 



If 5^("+iC,0) - ^|(K + */?) < 0, then we put ^ := 0, "+iCi := "Q, "+^Cii := '^C;; (no 
inelastic fiow occurs). Otherwise, ^ is computed using equation flHT]) . Substituting fl55]) o 
for in (IHTI) . we obtain two alternative forms of the incremental consistency condition: 



Dr'C,0-.= {^) '-^ ^ = 0. (89) 
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After the solution ^ is found, the values of ""'"^Ci, ""''^Cii are given by Ci("+^C,,^), 
Cii("+^C, ,^). Finally, we update s and Sd using equations fl82l) . 

Remark: Solving system flHT]) . fl83l) o with respect to C, with a given ^5^, it is possible to 
represent ^ as a function of '^'^^Ci, "''"^Cii, thus reducing the number of unknowns. On the 
other hand, for small 7] this approach will result in an ill-posed problem. 



3.4 Numerical implementation 



The Newton-Raphson method is used to compute (Ci("+^C, 0, Cii("+^C, O) from fl79|l. 
(IHOl) . To this end, equations (179|) . (IHOj) are linearized analytically using the coordinate- free 
tensor formalism proposed by Itskov (see [I6], |17j). 

Notice that the straightforward application of Newton's method to the solution of ( !88l) 
or ( |89l) is not trivial. Indeed, for r] = 0, m > 1, the convergence of the Newton method 
for (l88l) fails to be quadratic since the first derivative is zero at the root (see fig. [2]). At 
the same time, for i] > 0, the initial approximation ^'^"^ = can not be used to compute 
the solution of fl89l) . since the function -D(^) is not differentiable at zero (see fig. [2]). To 
overcome these difficulties, the first Newton iteration is performed using ( !88l) with initial 
approximation ^"^^^ = 0, and the subsequent iterations are performed using (l89l) . 

The derivative ^ — required by the Newton method, is calculated using the 

implicit differentiation of fl79|) . flHOj) with respect to ^. An alternative strategy is to 
solve fl89l) with the help of a derivative- free iteration scheme like Pegasus method [7], 
[20] . This approach is reasonable being combined with a fixed-point iteration for finding 
(Ci("+^C,0, Cii("+^C,0), such that no linearization of ([79]), ([80]) is required. 

We implement the coordinate-free tensor formalism to obtain an analytical expression for 
the consistent tangent operator ^n+i A ■ Using a special product A x B of two second-rank 
tensors and the composition A B of two fourth-rank tensors (see definitions (2.6), (2.10) 
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in [H]), it follows from (06]) that 

9"+^Ci("+^c) _ 9Ci("+^c,o 9Ci("+ic,o 5^("+^c) 

<9"+^C ~ ~ V di ) 9"+^C ■ 

The numerical computation of tensor exponential exp(B) is performed using Taylor power 
series expansion (!62l) . The derivative of tensor exponential is computed by (see [18j ) 

3exp(B)^g^g 

„=1 fc=o 

In general, this approach fails due to the roundoff errors, and more sophisticated tech- 
niques are required (see, for example, [30], [IE], [19], [25]). We do not use these advanced 
techniques in this paper, since in the present calculations the argument of the exponential 
function is bounded. Indeed, if ^ < 0.2 then it follows from (!85|) . (!86l) that 

ilBkll «2e<0.4, kG{i,ii}. (90) 

Therefore, the roundoff errors are negligible. Moreover, under condition (!90l) . the truncated 
power series only with few terms yield exact results up to machine precision. 



4 Numerical tests 



Now we analyse the accuracy of the integration methods presented in section 3. Toward 
this end, we simulate the material behaviour under strain controlled loading. The loading 
program in the time interval t G [0, 300] is defined by 



where 



with 



Fi := 1, F2 



F{t) = F'{t) or F{t) = F'(t), 

;i - t/100)Fi + (t/100)F2 if t G [0, 100] 

(2 - t/100)F2 + (t/100 - 1)F3 if t G (100, 200] 

(3 - t/100)F3 + (t/100 - 2)F4 if t G (200, 300] 



2 

1 

V2 




\ 



1 

V2J 



^10^ 



1 
1 







/i 

2 



V 









1 

V2J 



(91) 
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Remark. In this section we test the numerical schemes under a variety of loading con- 
ditions, in particular, under non-proportional loading. In this connection, the loading 
programm does not have to be mechanically plausible. 

The material parameters used in simulations are summarized in table [2l 



Table 2 

Material parameters 



k [MPa] 


fi [MPa] 


c [MPa] 


7 [MPa] 






73500 


28200 


3500 


460 






K [MPa] 


m [-] 7] 


[s-i] ko 


[Mpa] X 


[MPa^^] 


/?[-] 


270 


3.6 2 


10^ 1 


0.028 


5 



We put the following initial conditions on the internal variables 

Ci|t=o = 1, Cii|j=o = 1, s\t=o = 0, Sd\t=o = 0. (92) 



Only the uniform time stepping is used in this paper. The numerical solution obtained 
with extremely small time step {At = 0.01s) will be named the exact solution. 

The coordinates of Cauchy stress tensor T for loadings ( !9T|) ^ and ( !9TI) o are plotted respec- 
tively in figures [3] and m Note that det(F) = 1 if ( 19T1) ^ is used, and no hydrostatic stress 
occurs. On the other hand, relation ( l9Ti) o results in a large hydrostatic stress and a finite 
elastic bulk strain. 

The numerical simulation shows that both MEBM and EM have a similar error. Both 
methods produce slightly different results for At = 10 s when the inelastic increment ^ 
ranges up to about 17%. 



5 Characterization of the material model 



We investigate qualitatively the material response, predicted by the material model. The 
numerical computations simulate basic material testing experiments. Material parameters 
from table [2] and initial conditions (^21) are used in this section. 
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Fig. 3. Accuracy test with small elastic strains (use (f9T]) ^) 
5.1 Uniaxial testing 



1?00 



For uniaxial test we put 



^ 1 + e 0^ 



V 








a 
a 



-22 



33 



0. 



The unknown a is determined using (!93l) o. The technical stress a := ^^ii is plotted 
in figure Oa for various strain rates e. Here A, Aq denote the current and initial cross 



(93) 
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sections, respectively. Although the material response is stable, the stress reduction is 
observed after the peak load in uniaxial monotonic test. The reason is the reduction of 
the cross-section. The equilibrium curve can be reached both by relaxation and creep 
(figure Ob). In the simulation presented in figure Ob each relaxation period lasts for 
10 seconds. The creep time is 20 seconds. Therefore, the numerical experiment shows 
that it takes longer to reach equilibrium curve in the creep process than in the relaxation 
process. Finally, as indicated by the strain-controlled cyclic test (figure [51 c), the saturation 
is achieved after the isotropic hardening is accomplished. 
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0.15 
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0.04 0.08 0.12 0.16 0.2 

C 



Fig. 5. Uniaxial testing: monotonic loading (a), relaxation and creep (b), and cyclic loading (c). 



5.2 Torsion testing 



The cyclic torsion testing has much potential for providing information about the non- 



linear hardening phenomena (see 
tube we put 



[8]). To simulate the torsion of constrained thin- walled 



1 
ay 



33 



0. 



We consider a strain controlled torsion test with a given 
a is determined using (IMI) ^. Denote by a :— ^ 
stresses, respectively. 



T22 and T :- 



A 



(94) 



0.01/s. The unknown 
T12 the axial and shear 



The axial stress a is exactly zero in geometric linear theory. But, in the case of finite 
strains, so-called second order effects can appear, leading to nonzero axial stress. For 
instance, the Poynting effect (see [1]) is observed during the torsion of cylindrical samples 
made of aluminium alloy. This effect consists in axial compression of constrained samples 
or axial elongation of unconstrained samples (see \M^)- The axial and shear stresses are 
plotted in figure [6] for different forming increments. As may be seen from the figure, the 
Poynting effect is predicted by the material model. 

Next, as indicated by figure [6|, the maximal stresses are influenced by the forming incre- 
ment. For smaller forming increment the isotropic hardening is accomplished on the early 
stage of the forming process, thus leading to higher maximal stresses. On the other hand, 
if the kinematic hardening is not accomplished within one forming increment, a some- 
what different material response is possible. For instance, the maximal stresses under the 
cyclic loading can be smaller than the stresses under the monotonic loading. In order to 
demonstrate this effect, we perform the numerical simulation (see figure [7]) with modified 
hardening parameters: x := 0.0035 MPa"\ c := 1500 MPa, p := 10, 7 := 1800 MPa. A 
similar effect was reported in [27] for 20MoCr24 steel alloy. 
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Fig. 6. Torsion testing (material parameters from table [2]): shear stress (a), (b), and axial stress 
(c), (d). 

6 Discussion 



The classical material model of viscoplasticity is modified in a thermodynamically con- 
sistent manner to incorporate finite elastic and inelastic strains. The model takes rate- 
dependence (relaxation, creep) and hysteresis effects (nonlinear kinematic and isotropic 
hardening) into account. 

Although the material response is anisotropic, the symmetry of """"^Ci and ""''^Cii is a 
priori preserved by EBM, MEBM and EM. It is shown that no symmetrization procedure 
is necessary. Moreover, any symmetrization should leave the corresponding solutions un- 
changed. Both MEBM and EM have the advantage that the inelastic incompressibility 
constraint is exactly satisfied. 



Under special assumptions on the potential functions il^ei and V'kin it may be beneficial to 
optimize the solution procedure of system (1791) — (l86l) . thus reducing the computational 
effort. On the other hand, the most important properties of any stress algorithm are 
stability, accuracy, robustness and universality. The computational effort, required for the 
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Fig. 7. Torsion testing (modified material parameters): shear stress (a), (b), and axial stress (c), 
(d). 

evaluation of stresses and tangential operator, is negligible in comparison with the costs 
of solving the global linearized system of equations within the Newton-Raphson iterative 
procedure. 



The material model reproduces qualitatively the experimental results [15], [28] for alu- 
minium alloy processed by ECA-pressing. For more detailed modeling of kinematic hard- 
ening it is possible to introduce several Armstrong-Frederick terms, using series of mul- 
tiplicative decompositions of type ([3]). To complete the phenomenological description of 
the material, a proper parameter identification is required. 
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